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ABSTRACT 

Motivated by recent observations that show starless molecular cloud cores ex- 
hibit subsonic inward velocities, we revisit the collapse problem for polytropic gaseous 
spheres. In particular, we provide a generalized treatment of protostellar collapse in the 
spherical limit and find semi-analytic (self-similar) solutions, corresponding numerical 
solutions, and purely analytic calculations of the mass infall rates (the three approaches 
are in good agreement). This study focuses on collapse solutions that exhibit nonzero 
inward velocities at large radii, as observed in molecular cloud cores, and extends pre- 
vious work in four ways: (1) The initial conditions allow nonzero initial inward velocity. 
(2) The starting states can exceed the density of hydrostatic equilibrium so that the 
collapse itself can provide the observed inward motions. (3) We consider different equa- 
tions of state, especially those that are softer than isothermal. (4) We consider dynamic 
equations of state that are different from the effective equation of state that produces 
the initial density distribution. This work determines the infall rates over a wide range 
of parameter space, as characterized by four variables: the initial inward velocity foo, 
the overdensity A of the initial state, the index F of the static equation of state, and the 
index 7 of the dynamic equation of state. For the range of parameter space applicable 
to observed cores, the resulting infall rate is about a factor of two larger than found in 
previous theoretical studies (those with hydrostatic initial conditions and Voo = 0). 

Subject headings: hydrodynamics - stars: formation 
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1. INTRODUCTION 

Stars form through gravitational collapse and hence the dynamics of the collapse flow is funda- 
mental to understanding the star formation process. Historically, theoretical studies of gravitational 
collapse have considered initial conditions that correspond to hydrostatic equilibrium (or near equi- 
librium). In recent years, however, observations have shown that some molecular cloud cores display 
nonzero, inward velocities before they contain infrared sources. One standard interpretation of this 
finding is that the cores must begin their collapse with a head start, i.e., with a pre-existing inward 
flow. The goal of this paper is to reexamine the gravitational collapse problem for star formation 
with a focus on initial conditions that include initial inward flows, i.e., the conditions observed 
in star forming regions. Along the way, we generalize previous studies of protostellar collapse to 
include a wider set of equations of state. 

To isolate the effects of the inward flow, we return to the basic spherical collapse problem for 
polytropic spheres, which are used as idealized theoretical models (Shu 1977) of observed molecular 
cloud cores (e.g., Myers & Benson 1983, Benson & Myers 1987). This same treatment can be 
applied to the collapse of subcondensations (denoted as kernels - see Myers 1998) within larger cores 
that form stellar groups. In the simplest scenario, molecular cloud cores are thought to evolve into 
centrally concentrated configurations through the process of ambipolar diffusion (Mouschovias 1976; 
Shu 1983; Nakano 1984). Since recent observations suggest that core formation takes place more 
rapidly than predicted by standard models of ambipolar diffusion (e.g., Jijina, Myers, & Adams 
1999), several modifications to this picture have been suggested. These generalizations include 
ambipolar diffusion starting near the supercritical state (Ciolck & Basu 2001), turbulent fluctuations 
speeding up ambipolar diffusion rates (Fatuzzo &: Adams 2002; Zweibel 2002), core formation 
through turbulent cooling flows (Myers k. Lazarian 1998), and core (and star) formation through 
dynamic, turbulent fragmentation (Li, Klessen, & MacLow 2003; see also Li 1999). In the extreme 
limit of some efficient fragmentation scenarios, pre-collapse states (cores) never actually form at 
all. Our present approach - studying the collapse of cores that form with a pre-existing inward 
flow - is intermediate between the previous picture of star formation starting from hydrostatic 
equilibrium (Shu, Adams, h Lizano 1987) and the alternate scenario of a completely dynamic star 
formation process (MacLow & Klessen 2003). In any case, molecular cloud cores are observed to 
form relatively rapidly and to (often) exhibit subsonic inward motions in their starless state (Lee 
et al. 1999, 2001). These configurations thus specify the initial conditions for star formation and 
their subsequent collapse determines the manner in which stars form. 

The collapse of cloud cores can be studied using self-similar methods. Indeed, self-similar 
collapse solutions (Shu 1977) are an important cornerstone of the current theory of star formation 
(e.g., Shu et al. 1987). The original self-similar collapse calculations (Larson 1969ab; Penston 
1969ab; Shu 1977) considered an isothermal equation of state. Since then, many generalizations 
of the collapse have been made. The leading order effects of rotation have been studied, both 
for the inner pressure free region (Ulrich 1976; Cassen & Moosman 1981) and for the entire core 
(Terebey, Shu, & Cassen 1984). Different equations of state for the collapsing gas have been studied 
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(Cheng 1977, 1978; McLaughlin & Pudritz 1997). The leading order effects of magnetic fields have 
also been included (Galli & Shu 1993ab; Li & Shu 1996, 1997). More recently, the collapse of 
magnetized singular isothermal toroids has been studied (Allen, Shu, & Li 2003; Allen, Li, &: Shu 
2003); the latter study includes a calculation of collapse with an imposed initial velocity (see their 
Figure 7) , in the same spirit as this investigation. The effects of radiation pressure on the collapse 
flow have been studied in both the spherical limit (e.g., Wolfire & Cassinelli 1986, 1987) and the 
inner pressure free regime for rotating collapse (Jijina & Adams 1996). Additional mathematical 
generalizations of the solutions have also been considered (Hunter 1977; Whitworth & Summers 
1985). 

The self-similar collapse solutions explored in this paper are applicable for cores that exhibit a 
wide range of densities. Actual molecular cloud cores have a finite central density nc that extends 
over a given radial extent rc ~ as{2T:Giinc)~^/'^. The core must also have an outer radial scale 
Tout where the declining density profile of the core matches onto the background of the molecular 
cloud. If tm denotes the radius that encloses a stellar mass, then self-similar solutions apply over 
the ordering of radial scales rc <C tm ^ font- Previous numerical work (Foster & Chevalier 1993) 
indicates that the collapse of an isothermal core approaches the expected self-similar form when the 
core has Voutl'^'c > 20. For typical values M = 0.5 Mq and as = 0.2 km/s, this criterion implies a 
constraint on the central number density nc > 3 x 10'' cm~^. In the context of this paper, however, 
cores that have initial inward velocities can more readily approach the self-similar collapse forms 
and these constraints {rout/fc > 20; nc > 3 x 10^ cm~^) are less severe. 

Only a fraction of the observed starless cores have a sufficiently extended range of densities (or, 
equivalently, large enough ratios rout/fc) so that the collapse flows can be modeled to high accuracy 
with the self-similar forms derived here. Many observed cores have flat central regions (Tafalla et 
al. 2004), and tend to approach the form p ~ r^^ beyond a few thousand AU (Ward-Thompson, 
Motte, & Andre 1999). In a recent survey (Bacmann et al. 2000), 3 of the 24 observed cores also 
show evidence for a well defined edge (corresponding to the scale Tout) so that these cores have 
small ratios Vout/rc- However, mapping results in Taurus are consistent with a self-similar collapse 
picture and the envelopes in Bok globules are also qualitatively consistent (Motte & Andre 2001). 
For cores with constant density centers, the infall rate will display more time variability than in the 
self-similar solutions: The infall rate will be very small at first, but quickly grows to become larger 
than that of the self-similar solution; once the constant density portion of the core has collapsed, 
the infall rate decreases toward the self-similar value (e.g., Foster & Chevalier 1993; Myers 2004). 
Furthermore, even if actual cores only approach the self-similar solutions asymptotically in time, 
this idealized formulation of the collapse problem allows us to isolate the effects of various physical 
inputs (especially nonzero initial velocities). 

The overall goal of this paper is thus to obtain a greater understanding of the collapse of 
molecular cloud cores. This work includes four generalizations of previous studies of protostellar 
collapse, but the most important astronomical issue is to include the effects of pre-existing inward 
fiows, which have now been observed (as outlined above). We can include an initial inward velocity 
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as part of the initial conditions (i.e., we let v ^ at large radii). In a similar vein, we can consider 
initial states that are overdense and thus exceed the density required for hydrostatic equilibrium. 
These two complications are related in that both can explain the observed inward flows and both 
lead to larger mass infall rates, but they are conceptually distinct. In the first case, the inward 
flows are part of the initial conditions and must be produced by whatever physical processes form 
the core. In the case of overdense initial states, however, the inward flows are part of the collapse. 
Specifically, in the latter case, the core begins with zero initial velocity throughout its structure, 
but the overdense core soon develops small, but nonzero, inward speeds at large radii (as observed). 

In addition, we consider equations of state that are softer than isothermal. This generalization 
is motivated by observations of non-thermal linewidths Av in molecular cloud cores, where one 
finds correlations with density of the form (Af)^ P^^ ■, with /3 > (e.g., Larson 1985, Jijina et 
al. 1999). If one interprets the linewidth as the effective transport speed in the medium, then 
the corresponding effective equation of state has the form P ^ , where V = 1 — (5 and hence 
< r < 1. In the limit T ^ P ^ logp; this limiting case of a sequence of possible equations of 
state is often called a "logatrope" (beginning with Lizano k. Shu 1989). Here we consider a range 
of effective equations of state, from the isothermal limit F = 1 down to the logatropic limit F ^ 
(collapse solutions in the logatropic limit have been studied previously by McLaughlin &: Pudritz 
1997). Notice that the general equations we work with below reduce to the logatropic case in the 
limit r ^ (but care must be exercised in taking the limit). 

As another generalization, we consider the case in which the dynamic equation of state for 
the collapsing gas is different than the effective (static) equation of state that produces the initial 
equilibrium configuration. Here, the static equation of state refers to the pressure law that enforces 
the initial (pre-collapse) configuration for the core. Note that we consider only barotropic equations 
of state throughout this paper. The dynamic equation of state (as set by 7) refers to the pressure 
law that describes how the thermodynamic variables of the gas change as the material is compressed 
during collapse; this process is governed by the entropy evolution equation (see §2.1). Equations of 
state that are softer than isothermal are motivated by the observed linewidth dependence on density, 
and the inferred density profiles, in molecular cloud regions. These considerations apply primarily 
to the static equation of state as defined here (for applications of such soft, static equations of state, 
see, e.g., Curry & McKee 2000, Curry & Stabler 2001, Spaans & Silk 2000). The dynamic equation 
of state, as set by the index 7 that appears in equation (5), need not be the same as the static 
equation of state. In other words, the physics that determines the density profiles of the pre-collapse 
states can be different from the physics that governs the thermodynamics of the collapse flow. We 
thus allow 7 7^ r. This possibility is also motivated by theoretical MHD simulations (Vazquez- 
Semadeni, Canto, & Lizano 1998) which show that the turbulent velocity dispersion increases with 
mean density as the gas collapses, in contradiction to the expected behavior for a soft equation 
of state; as a result, logatropic and similarly soft equations of state may provide an inadequate 
description of the dynamical processes occurring within a cloud. 

The mass infall rate M is one of the most important physical quantities in the star formation 
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problem and this paper focuses on calculating M for a range of conditions. Collapse flows are often 
self-similar and have no characteristic mass scale. Instead, the collapse flow feeds material onto the 
central star /disk system at a well-defined mass infall rate M. In these flows, the infalling material 
always approaches free-fall conditions on the inside (in the limit r ^ 0) and the reduced mass mo 
determines the size of the infall rate (as shown below). The infall rate M determines, in part, the 
total system luminosity and the total column density of the infalling envelope. These quantities, 
in turn, largely account for the spectral appearance of protostellar objects (e.g., Adams, Lada, & 
Shu 1987; Adams 1990). For protostellar objects, most of the luminosity is derived from material 
falling through the gravitational potential well of the star. Although the circumstcllar disk stores 
some of the energy in rotational motion, the system luminosity is (usually) a substantial fraction 
of the total available luminosity 

L.^'^. (1) 

-ft* 

where M is the total mass of the system, M is the mass infall rate, and i?^ is the stellar radius. 
The stellar radius, which helps determine the depth of the potential well, is itself a function of the 
mass infall rate (Stabler, Shu, & Taam 1980). 

This paper is organized as follows. In §2, we present a general formulation of the collapse 
problem, including both the similarity transformation for self-similar solutions and a complementary 
numerical treatment. We present general analytic results in §3, including a discussion of how the 
"outer" spherically symmetric solutions of this paper match onto collapse solutions in the inner 
regime where the effects of rotation are important. In §4, we consider collapse solutions for which 
the static equation of state is isothermal (so that the initial state p ^ r~^), but generalize previous 
work to include nonzero initial flow speeds [v^ > 0), overdense initial configurations (A > 1), and 
dynamic equations of state that are not isothermal (T = 1, 7 7^ 1). In §5 we consider collapse 
solutions from non- isothermal initial conditions, i.e., starting states for which the static F 7^ 1. In 
§6, we conclude with a summary and discussion of our results. 



2. FORMULATION OF THE COLLAPSE PROBLEM 

In this section we formulate the collapse problem for the general case of a collapsing spherical 
cloud of gas. The collapse solution is obtained by self-similar methods (§2.2) and by a numerical 
approach (§2.3). In both cases, we allow the initial states to be out of exact hydrostatic equilibrium 
by having initial inward velocities and/or by being overdense. We also allow the dynamic equation of 
state to be different from the static equation of state that determines the initial density distribution. 
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2.1. Basic Governing Equations 

The basic equations governing the fluid are given below. We want to consider the spherical 
collapse problem, so that the equations of motion reduce to a simpler form. Conservation of mass 
can be expressed in terms of mass shells 

dM dM ^ , ^ 2 

— + u— = and ^ = Wp, (2) 

where M{r,t) is the total mass within a given radius r at a time t. The above equations are 
equivalent to the usual equation of continuity, 

dp I d , r, ^ 

where p is the density and u is the radial component of the velocity. The force equation can be 
written in the form 

du 

dt dr p dr r 
where P is the pressure. Finally, the conservation of entropy equation can be written in the form 



+ ^^ = -T^-^. (4) 



(| + ^|:)[iog(P/p^)] =0, (5) 

where 7 is the index of the dynamic equation of state. Adopting equation (5) is not the same as 
assuming an equation of state of the form P = Kp"^. Equation (5) describes how a given parcel of gas 
changes its thermodynamic variables along a streamline, whereas the equation of state {P = Kp'^) 
implies a global constraint on those variables. 



2.2. The Similarity Transformation 

In this section, we find a similarity transformation for the collapse problem formulated above. 
Mathematically, this problem is represented by coupled partial difl'erential equations in the variables 
time t and radial position r. Using the similarity transformation, wc can reduce this problem to 
a set of ordinary difl'erential equations in a new similarity variable x which we deflne below. In 
particular, we look for a similarity transformation of the general form 

X = Afr , p = Bt^a{x) , M = Cfm{x) , 

u = Dt'^v{x), and P = Et^p{x) . (6) 

Here, both the coefficients [A, B, C, D, E) and the indices (a, h, c, d, e) are constants. The reduced 
fluid fields (a, m, v, p) are dimensionless functions of the (single) dimensionless similarity variable 
X. The time benchmark t = corresponds to the instant of the onset of collapse when the mass of 
the central "object" is zero, i.e., M{0,t) = 0. 
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The general similarity transformation calculation leads to four equations to specify the five 
indices o, b, c, d, e. We leave the constant a arbitrary for the moment and write the rest of the 
variables in terms of its value, i.e., 

a = a, b = —2 , c=— (3a + 2), 

d = -(a + l), and e = -2(a + 2) . (7) 
Similarly, for the coefficients we obtain 

A = A, B= (47rG)-^ , C = {A^G)-^ , 

D = A-^, and E = {A-kGA^Y^ . (8) 

Keep in mind that G is the gravitational constant. We thus obtain reduced equations of motion in 
the form 

din 

(ax + v)— = (3a + 2)m , (9) 

dx 

dm 9 , , 

— =x^a, (10) 
dx 

.dv 1 dp m , ^ , ^ 

{ax + v)— + -^ = ^ + {a + l)v, (11) 

dx a dx x^ 

[ax + v)4- logb/«^] = 2(2 + a - 7) , (12) 
dx 

, , 1 da dv ^.^ , ^ „, 

(ax + v)-— + — = 2{l- V x) . (13) 
a dx dx 

This similarity transformation is not unique - one can always rescale the coefficients {A, B, C, D, E} 
by a set of dimensionless numbers and obtain new equations of motion with different numerical 
coefficients. 

Notice that wc can immediately combine the first two equations of motion to obtain an ex- 
pression for the reduced mass m{x), i.e., 

(ax + v) 2 

m = -. -X a . (14) 

(3a + 2) ^ ' 

Next we define a constant q according to 

If we then take q times equation (9) and subtract it from equation (12), we can integrate the 
resulting equation to obtain an expression for the reduced pressure p{x), i.e., 

p = Coa^+^ {x^ + vx'^/ay = Ci m« , (16) 

where Cq and Ci are constants which depend on the initial density distribution of the core. Given 
these solutions for the reduced pressure p(x) and reduced mass m{x), equations (11) and (13) are 
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the relevant equations of motion to determine the remaining unknown functions a{x) and v{x). 
Notice that in the hmit x oo, the sohition (16) for the pressure approaches the form p C2cf 
(obtained by integrating the equations of motion). In the opposite limit x ^ 0, p ^ C^a^ . In other 
words, this solution for the pressure illustrates how the equation of state smoothly transforms from 
static index V (at early times and/or large radii) to the dynamic index 7 (at late times and/or 
small radii). 

In this paper, we consider only positive values of time, since we take the zero point of the 
time variable to be the onset of collapse. It is mathematically possible to consider solutions for 
negative times, t < 0, which would correspond to negative values of x. Solutions that span the 
entire available range — 00 < x < 00 are known as "complete" solutions (first obtained by Hunter 
1977; see also Whitworth & Summers 1985). Note that the self-similar solutions of this paper for 
the protostellar collapse phase {t > 0) cannot (in general) be extended to the pre-stellar phase 
{t < 0), as they would likely encounter a critical point and become singular. In other words, 
dynamical evolution (using only collapse physics) in the t < regime would not (in general) lead 
to the initial conditions used here. However, it is rather unlikely that molecular clouds will evolve 
toward their centrally condensed initial configurations in a self-similar manner subject to (only) the 
physics included in these equations of motion. Before the onset of collapse (for t < 0), molecular 
clouds may evolve through the processes of ambipolar diffusion (at least for small mass scales), 
shocks, turbulent dissipation, cooling flows, condensation instabilities, and cloud-cloud collisions. 
In addition, the cloud will most likely initiate collapse before a completely self-similar equilibrium 
state has been attained; the collapse will only become self-similar asymptotically in time (i.e., the 
self-similar collapse solutions of this paper are intermediate asymptotic solutions to the realistic 
problem of the collapse of a finite cloud with finite central density). In any case, we limit this 
discussion to solutions with < a; < 00, sometimes called "semi-complete" solutions. 



2.3. Numerical Treatment 

As a complement to the self-similar formulation, we also perform numerical simulations of the 
collapse for the case where F = 7 = 1. For the isothermal equations of state considered here, 
P = a^p, the force equation can be written as 

^(H + -2Q-^irpu) = -a,- - p^ . (17) 

This form of the force equation is solved numerically along with equation (3) using the second- 
order-accurate scheme outlined in Boss Sz Myhill (1992). This scheme invokes a variation on the 

standard predictor-corrector technique, and allows for nonuniform grid spacing. Our grid is defined 
by 160 logrithmically spaced radii (r^) spanning four orders of magnitude from the inner to the 
outer boundaries. The cell interfaces (^^^1/2) ^^'^ then set at the midway points of these radii. The 
mass terms at each radius (Mj) are advanced in time (denoted by the superscript n) using the 
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scheme 



where the flux terms are defined as 



n n 



at ' \dt^ 

and the right hand terms in the above equations are defined in Boss & Myhill (1992). 



(18) 



(19) 



3. GENERAL ANALYTIC RESULTS 

Given the formulation of the collapse problem, we now find several analytic results. These 
results are general in that they apply to the entire class of solutions considered in this paper. 



3.1. Relationship Between the Similarity Transformation 
and the Static Equation of State 

First, we show that the initial density profile of the core determines the type of similarity 
transformation which describes the subsequent collapse of the core. To start, we assume that the 
initial equilibrium configuration of the system can be described by a "static" equation of state of 
the form 

P = Kp^ . (20) 

Since the similarity transformation determines the manner in which both the pressure P and the 
density p must scale with time, this new equation of state imposes an additional constraint on the 
system. In particular, the similarity transformation requires that the scaling exponent a (which 
has been left arbitrary thus far), must have the value 

a = r-2, (21) 

where the similarity variable is defined hj x = At^r. We also obtain a required value for the 
corresponding numerical coefficient A, i.e., 

A = K-V2(4;rG)(r-i)/2 . (22) 

With this choice for A, the reduced static equation of state has the form p = . As mentioned 
before, however, this coefficient A can be rescaled by a dimensionless constant without changing 
the physics (as long as that constant is propagated throughout the set of equations of motion). 

Notice that for the special case in which the static and dynamic equations of state are the 
same (i.e., F = 7), the index q = (see equation [15]). In this case, Cq = 1 = C\ and the equation of 
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motion which defines the pressure has the obvious solution p = oF . In the general case, the index 
q is given by 

«=^^- (23) 

^ 3r-4 ^ ' 



3.2. Hydrostatic Equilibrium and Overdensity 



One useful reference state is that in which the core is in hydrostatic equilibrium at t = when 
the collapse begins. In the limit of large values of we thus expect the core to have no velocity so 
that V = In this limit, x 3> the equations of motion have solutions of the form 



'm{x) = 



2-r 
4-3r 



2r(4 - 3r) 
(2 - vy 



i/(2-r) 



,(4-3r)/(2-r) 



a[x] 



2r(4 - 3r) 
(2 - r)2 



i/(2-r) 



X 



-2/(2-r) 



(24) 



(25) 



The density profile of the hydrostatic equilibrium configuration is thus a power-law with index 
Jl= —2/(2 — r). Notice that the static equation of state alone determines the initial configuration 
(as expected). Notice also that solutions of this form make sense only for static equations of state 
with r < 4/3. This latter result is, of course, well known from stellar structure theory (e.g., 
Chandrasekhar 1939). 

In this paper, we allow for the initial states to have densities (and mass profiles) that are 
heavier than that appropriate for hydrostatic equilibrium. For general (static) equations of state, 
we let the parameter A denote the overdensity, so that the initial mass and density profiles are 
given by the above equations with a leading factor of A > 1. For an isothermal (static) equation 
of state, however, we follow previous authors (e.g., Shu 1977) and write p(r) = Aal/2'KGr'^ so that 
the starting profiles have the known forms a{x) = A/x^ and m{x) = Ax. In this case, A = 2 
corresponds to hydrostatic equilibrium and A > 2 represents overdense initial states. Note that 
K = A/2 iox isothermal initial configurations. 

It is useful to consider the density and mass profiles in terms of physical quantities. Including 
the constants that carry dimensions and the overdensity parameter, we can write the profiles in the 
form 



M{r) = 47rA- 



3r 



kF 4 - 3F 



i/(2-r) 



and 



p(r) = A 



27rG (2 - F)2 

i/(2-r) 



,(4-3r)/(2-r) 



kT 4-3F 
2ttG (2 - r)2 



-2/(2-r) 



(26) 



(27) 
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3.3. Location of the Expansion Wave 

We can derive a general expression for the location of the expansion wave. If the initial state of 
the cloud core is in hydrostatic equilibrium, then the boundary between the inner collapsing region 
and the outer static region will be a critical point in the flow. The mathematical manifestation of 
this critical point is that the matrix of coefficients in the system of differential equations (9 - 13) 
must vanish (at the critical point). This condition can be written in the general form 

j = a^-^i2-rfxl, (28) 

where we have used the fact that v = and p = at the head of the expansion wave xh- Using 
the solutions for the starting density profile derived above, we can eliminate the reduced density a 
and find the following expression for the location xh of the head of the expansion wave 

XH = 7(1-^/2) (2 - r)-i [2r(4 - 3r)]-(^-^)/^ . (29) 

Notice that this expression is valid for initial states that are in hydrostatic equilibrium (A = 1) and 
does not include the terms appropriate for nonzero initial velocities. In our analytic derivations, 
we use this expression as a useful benchmark for estimating infall rates. For future reference, it is 
also useful to have the (total) mass mn enclosed within the expansion wave. Using the starting 
density profile and the expression for xh, we find 

niH = 2r7(2-3r/2) [2r(4 - 3r)]-'^^^-^^'^^^{2 - F)-^ . (30) 

3.4. Time Dependence of the Mass Infall Rate 

The static equation of state also determines the time dependence of the mass infall rate M. 
The infall rate determines the time scales for the star formation processes. In addition, since these 
cores have no well defined mass scales, the mass infall rate is the physical quantity of importance 
for star formation. Using the above results, we find that the mass infall rate has the form 

M = mo k'/^ (4 - 3r) (47rG)3(i-r)/2 q-' i^d-r) ^ ^3^) 

where mo is a constant which represents the reduced mass m{x) in the inner limit x 0. Since 
the infalling material always approaches ballistic conditions in this inner limit, the magnitude of 
the reduced mass mo determines the infall rate. Notice that, except for this parameter niQ, which 
is expected to be roughly of order unity, the mass infall rate is completely specified by our scaling 
transformation. The remainder of this paper is (mostly) devoted to the determination of mo for 
the cases of interest. 

For the case of an isothermal initial state, F = 1, and the infall rate M is constant in time. 
Notice that for equations of state softer than the isothermal case (F < 1), the mass infall rate 
increases with time. On the other hand, for harder equations of state (F > 1), the mass infall rate 
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decreases with time. Since we expect most star forming regions to have initial core configurations 
corresponding to soft equations of state, the general trend is for mass infall rates to be increasing 
functions of time. 



3.5. Matching onto Inner Solutions: Effects of Rotation 

In this section, we determine the most important effects of rotation on the collapse solutions. 
For the cases of interest, the inner limit of the spherical collapse flow approaches a ballistic (pressure- 
free) form. However, the spherical approximation must breaJc down in the inner region of the flow 
where conservation of angular momentum plays an important role and where a circumstellar disk 
forms. We thus consider the inner limit of our spherical solutions as the outer limit of the rotating 
non-spherical solutions which we calculate below. This exercise in intermediate asymptotic analysis 
remains valid under the following ordering of (radial) size scales: 

i?* < i?C < r-if . (32) 

In this ordering constraint, the scale ii* is the radius of the forming star and defines the inner 
boundary of the collapse flow. The scale Rc is the centrifugal radius (defined more precisely 
below) which roughly divides the spherical outer region of the flow from the highly nonspherical 
inner region. Finally, the scale rn is the head of the expansion wave and (again roughly) divides 
the outer core from the collapsing inner core. This latter division is clean only for hydrostatic 
initial states. 

Here, we present the form of the density profile resulting from the collapse of slowly rotating 
cloud cores. This calculation generalizes the isothermal case studied previously (Tcrcbcy, Shu, & 
Cassen 1984), where the collapse solution in the outer portion of the core matches smoothly onto 
an "inner solution" that can be described analytically (Cassen k, Moosman 1981 and Ulrich 1976). 
In this inner regime, parcels of gas spiral inward on nearly ballistic trajectories with nearly zero 
total energy (and conserved specific angular momentum). The resulting orbits are parabolic and 
are described by the equation 

/^o - ^ C 

where /xq is the cosine of the angle of the asymptotically radial streamline (i.e., the parabolic 
fluid trajectory that is currently passing through the position given by C, and ji = cos 9 initially 
made the angle with respect to the rotation axis). Although equation (33) is cubic in /xq and 
has three roots, only one solution has physical signiflcance. The quantity ( is defined by 

C = = — (34:) 

^ - GMr r ' ^ ' 

where joo is the specific angular momentum of parcels of gas currently arriving at the origin along 
the equatorial plane. We have followed previous authors in assuming an initial state that is rotating 
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at a constant rotation rate J7, so that the quantity joo is given by 



3c 



(35) 



where is the starting radius of the material that is arriving at the origin at a given time. This 
radius can be determined by inverting the initial mass distribution M(r), i.e., 



roo(M) 



M 
4ttA 



(2-r)/(4-3r) 



2ttG 



(2 - r)^ (4 - 3r) 



i-r 



i/(4-3r) 



(36) 



Putting all of these results together, we can find the centrifugal radius for any given static 
equation of state. In general form, the centrifugal radius can be written 



Rc = 



^-4(2-r)^4-r^3(r-i)^2Kr)-^(4 - 3Ty-''^i2 - Vf'^-^inf'-'] ^'^^ , (37) 



where the parameter A > 1 is the factor by which the initial states are denser than that required 
for hydrostatic equilibrium. The result for three specific cases are given below. For isothermal 
equations of state, F = 1, we obtain the familiar result 

Rc = 5 , (38) 

in the logatropic limit 

^ = (30) 

where Pq is the coefficient in the equation of state (which has units of pressure in this limit). 

Using equation (33) in conjunction with conservation of angular momentum and conservation 
of energy one can determine the velocity field (Cassen Sz. Moosman 1981). The density distribution 
of the infalling material can be obtained by applying conservation of mass along a streamtube 
(Terebey, Shu, & Cassen 1984; Chevalier 1983), i.e., 

p{r, e) Vr sin ed9d(f) = -^ sin 9o d9o dcf). (40) 

Combining the above equations allows us to write the density profile in the form 

-1/2 r 1-1 



n(r 6) = - 1 + -^ 



l + 2CP2(/^o) 



(41) 



where ^2(^0) = (S/Xq ~ is the Legendre polynomial of order two. Notice that this form for the 
density distribution is general in that it applies to all of the collapse calculations considered in this 
paper. The properties of the collapsing core determine the form of the functions M and ( = Rc / r 
which appear in the general form (40). 



-14- 



4. COLLAPSE SOLUTIONS FROM ISOTHERMAL INITIAL STATES 

The solution for the cohapse of a singular isothermal sphere has been found previously (Shu 
1977) and has been used in numerous applications of star formation theory. Most applications 
use the "expansion-wave collapse solution" which starts from a hydrostatic initial state and has 
zero initial velocity. To account for the observed inward flows in starless cores, one can generalize 
the solution to include overdense initial states (as considered in Shu 1977) and/or nonzero starting 
velocities. In this section, we include both of these complications and show the relationship between 
them. We also generalize the collapse solution to include the effects of a dynamic equation of state 
which differs from isothermal, i.e., we set F = 1 but allow 7 7^ 1. For these solutions, the parameter 
K that appears in the static equation of state has units of speed squared and can be written in 
terms of the isothermal sound speed, i.e., 1^ = 0^. 



4.1. Self-similar Solutions 

For the case of an isothermal initial state, with static F = 1, the parameter q = 2(7 — 1) and 
the reduced equations of motion take the form 

, ,dv 1 dp m , 
{v-x)— + --f = -—, 42 
ax a ax x^ 

, . 1 da dv , , , 

where the reduced pressure p{x) and mass m(x) are given by 

p = Ca"' m'^^'^~^^ and m = x^a{x — v) . (44) 

In this paper, we explicitly consider the starting states to exhibit two generalizations from 
exact hydrostatic equilibrium. The cores can be overdense, i.e., with densities larger than can be 
pressure supported so that A > 2. We also consider the possibility of nonzero initial flows as 
characterized by the initial speed v^o- As a result, we take the starting states to have the form 



A 

a = ^ 

and 



l-^ + 0{x-') , (45) 



2x2 

A-2 (^ _ 2) [(5 - 27)/3 - A/6] - 2vl , , 4 



X x^ x^ 



+ 0(x-^) . (46) 



The density profile has the same form as considered previously for overdense states (Shu 1977), 
whereas the velocity field picks up additional terms due to 7^ and 7 7^ 1. The limit of 
hydrostatic equilibrium corresponds \o A^2 and Voo — Oj so that v — >^ and a^2jx^. 

For isothermal starting states (static F = 1), we are left with a three parameter family of 
collapse solutions, as specified by the density parameter A^ the inward velocity scale and the 
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index 7 of the dynamic equation of state. For each choice {A,Vcx,,^), we get a different similarity 
solution, and, in particular, a value of the reduced mass at the origin ttt-q (which specifies the mass 
infall rate). The standard infall-collapse solution (Shu 1977) has rriQ = 0.975 and corresponds to 
a particular point in this parameter space, namely (2,0,1); that same paper also presents solutions 
for points {A, 0, 1), although they are not widely used in subsequent work. 

The first results of this paper are summarized by Figures 1 and 2. In Figure 1, we plot the 
reduced mass mo as a function of the overdensity parameter A for various values of the index 7 
appearing in the dynamic equation of state. Figure 1 shows that the dynamic equation of state has 
only a modest effect on the mass infall rate, whereas variations in the overdensity A are much more 
important. Similarly, Figure 2 shows the variation of the reduced mass mo as the initial inward 
velocity Voo increases. In this case, the dynamic equation of state affects the niQ values more than 
in the case of overdense states, but variation of the inward speed has greater influence. 

4.2. Numerical Collapse Solutions 

The numerical treatment of the collapse problem can be used to verify the semi-analytic re- 
sults and to study alternate collapse scenarios that are not fully self-similar. For this numerical 
exploration, we specialize to the case F = 1 = 7. We consider collapse solutions for various values 
of the overdensity parameter A, where A = 2 is the value appropriate for hydrostatic equilibrium 
(where we adopt the notation of Shu 1977 for isothermal starting states). In addition, we allow for 
a non-zero infall velocity Uin in the initial state. As discussed previously, our self-similar treatment 
considers nonzero velocities in the outer region with the specific profile given by equation (46) for 
isothermal initial conditions. We note that current observations measure the inward velocities at 
a limited range of local radii, so a variety of different velocity fields (inward speed as a function of 
radius) remain possible. Here, we numerically determine the infall rates for a collection of isother- 
mal spheres with nonzero outer velocities that are constant (with radius). The speed nj„ and the 
overdensity parameter A are thus the two parameters that can be varied. 

Figure 3 shows one representative model with A = 2.0 (a density profile in hydrostatic equi- 
librium) and initial inward speed Uin = 0,3 /2 (a typical observed value). The top panel of Figure 
3 shows the density profile at the initial time (dotted curve) and four subsequent times. Except 
for the overall normalization constant, this solution is nearly identical to that found in the original 
self-similar collapse solutions (compare with Figure 3 of Shu 1977). The second panel shows the 
velocity profile. The inner regime approaches the ballistic (free-fall) form, \v\ ~ r~^/^, as expected. 
On the outside (large radii), the inward velocity attaches smoothly onto the imposed boundary 
condition v = constant. The numerically calculated infall rate quickly approaches the constant 
value predicted by the similarity solution and stays fixed at that value. 

Figure 4 illustrates that the flow is indeed self-similar (as expected). Here, the collapse solutions 
for both the density (top panel) and velocity field (bottom panel) at four different times are rescaled 
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according to the similarity transformation of §2 and then plotted together. The result is one smooth 
function for both density and velocity, as expected for self-similar flow fields. 

The resulting values of iuq are shown in Figure 5 for a set of models with varying overdensities 
(specified by A) and varying infall speeds in the range Voo = Uin/O's = - 1.0. Figure 5 shows 
that the mass infall rates (specified through mo) increase with increasing overdcnsity A and with 
increasing initial inward speed v^o- Furthermore, the functional dependence of mo on both A and 
Voo is approximately linear so that increasing A and increasing Voo produce similar effects on the 
collapse solutions. A more detailed comparison of the effects of A and Voo on the mass infall rate 
emerges from the analytical estimates derived in the following subsection. 



4.3. Analytic Estimates 

In this subsection, we derive analytic estimates for mo for the case of initial configurations 
with static F = 1. We consider, separately, both the case of overdense initial states with A> 2 and 
initial states with nonzero inward velocities Voo at the start. To start, for simplicity, we set dynamic 
7 = 1. These estimates, by necessity, require approximations, but in the end they provide both 
analytic understanding of the collapse physics and accurate analytic formulas for the integration 
constants mo and hence the infall rates. 

We first consider the case of overdense initial states. For a given fluid layer at initial radius 
ro, the inward velocity is zero at t = (where x — ^ oo), but an inward flow soon develops. Due to 
the self-similar nature of the problem, the value of the radius ro does not matter. To estimate mo, 
we assume here that the fluid layer migrates slowly inward as a nonzero inward velocity develops. 
As the layer moves inward, the expansion wave moves outward. When the expansion wave passes 
the fluid layer on its way out, we assume that the layer subsequently falls inward in a pressure- free 
manner. For a starting radius ro, the fluid layer has a radial location at later times given by 

r^ = rl-{A-2)alt\ (47) 

where Qg is the sound speed of the initial state. The expansion wave has radial location rn = cist, 
so the expansion wave crosses the fluid layer at time tc given by 

te = (ro/a,)(^-l)-V2. (48) 
At this time, the layer resides at radius rc given by 

rc = rH = astc = ro{A-l)-'^/^ . (49) 

From this location at rc, the layer collapses in a pressure-free manner. The layer is already moving 
inward with speed 

Uino = ^^^^—^t, = aM - 2) , (50) 
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where we have kept only the leading order term in the expansion. With this starting speed, the 
time required for the fluid layer to fall to the origin is given by 

tf = r'J'[2GM{ror'/'m, (51) 

where we have defined 

fin) ^ (1 - r?)-3/2{sin-\l - r/)V2 _ (^ _ ^2)i/2| ^ (53) 

where the parameter rj = uf^^.rc/[2GM(ro)] measures the size of the initial velocity. The effect of 
Uinc is thus determined by the dimensionless function /(r/). Here, the enclosed mass is M(ro), i.e., 
all of the mass within the starting radius ro, even though the layer now resides at the smaller radius 
Tc. As result, the relevant value of 77 is given by 

{2A){A-iy/'- ^^^^ 

In the limit of hydrostatic equilibrium, A ^ 2 and rj ^ 0, so that we recover the classical result 
/(r/) — 7r/2. Since the inward speed cannot be larger than the free-fall value, this formula is only 
valid for A < Amax ~ 9.35 (although we are interested in more modest values of A). The infall 
time can be written 

The total time for the fluid layer to fall is then given by the sum It = tc + tf, i.e., 

rp/as f{ri)ro/as , . 

^ {A - 1)V2 ^ (A- 1)3/4(2^)1/2 • ^^^^ 

By the time tr that the fluid layer reaches the origin, the expansion wave has traveled out to 
rn = dstr- After some algebra, the ratio of reduced masses can be written in the form 

rriH ^ J 

Since mn = A, this estimate for mo is now completely specified. Notice that for sufficiently large 
value of A (namely, .4^ 2.85), the ratio mo/iTiH can exceed unity. This occurs because the inward 
velocities outside the nominal location of the expansion wave move material inward faster than the 
outward-moving expansion wave can enclose more material. 

Figure 6 shows a comparison of our self-similar formulation, the analytic estimate derived 
above, and results from our numerical treatment. For this comparison, all initial states are taken 
to be isothermal (so that static F = 1) and the collapse flow is assumed to be isothermal (so that 
dynamic 7 = 1). The resulting values of the dimensionless mass mo, which specifies the mass infall 
rate, are shown as a function of the overdensity parameter A. Notice that all of the curves are in 
good agreement. The self-similar results (solid curve) have been determined previously (see Table I 
of Shu 1977), but such overdense solutions have been largely ignored. The numerical results (dotted 
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curve) demonstrate that the self-similar formulation produces the correct infall rates for cores that 
start with the same initial conditions (here, density profiles that are sufficiently centrally condensed 
- sec also §4.4). In addition, the analytic estimates are in excellent agreement with both of the 
more rigorously determined values. The analytic treatment makes two approximations, which tend 
to cancel out: We first assume that once the expansion wave overtakes a fluid layer, the material 
immediately becomes ballistic, and hence this assumption acts to overestimate the infall rate. We 
also assume that the initial speed in given by equation (50) , which keeps only the leading term in a 
series expansion, and hence tends to underestimate the infall rate. As a result, the dependence of 
niQ on the overdensity (A) is well described by equation (56). In spite of its derivation, the result 
(eq. [56]) provides an accurate and user-friendly formula to specify infall rates. 

We now present an analogous calculation for the case of initial inward velocities (ttj„ 7^ 0). To 
isolate the effects of the velocity, we fix the density coefficient to be that appropriate for hydrostatic 
equilibrium so that A = 2. For a fluid layer at starting radius ro, the inward speed is nearly constant 
at Uin, and the radial location as a function of time is given by 

r{t) = ro - Uint . (57) 

Using the standard expression for the expansion wave, which is only an approximation in the present 
context (see eq. [28]), we can find the time tc at which the expansion wave crosses the fiuid layer, 

tc = ro{uin + as)~^ , (58) 

and the corresponding location Tc of the layer, 

rc = rQas{uin + as)~^ . (59) 

As before, we assume that the fluid layer collapses subsequently in a pressure-free manner. The 
additional time tf (beyond tc) required for the fiuid layer to reach the center is still given by 
equation (54), with the intermediate crossover radius rc given by equation (59) and the parameter 
rj given by 

Keep in mind that v^o = Uin/ag. Thus far, we have assumed that the fluid layer moves inward 
at constant speed Uj„ until it crosses the expansion wave and then falls inward with no pressure 
forces. In actuality, the layer will accelerate a small amount before crossing the expansion wave, 
but will be held back (again, a small amount) by pressure forces; these two corrections thus tend to 
cancel. However, we can account for this complication by introducing a dimensionless parameter b 
through the ansatz = hvoo- After some algebra, the ratio of reduced masses takes the form 



^={[^i4/[„(5^)i[ ' (61) 



The mass contained within the expansion wave mn = 2 for these conflgurations, so the mass mo is 
completely specifled. 
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Figure 7 shows how this analytic formula compares with the results of the self-similar formu- 
lation and the numerical treatment of the collapse problem. All of the initial states are taken to 
have ^ = 2, as well as F = 1 =7. The resulting values of the reduced mass niQ are shown as a 
function of the starting inward speed v^^ = Uin/ag- The solid curve shows the results of the self- 
similar calculation and the dashed curve shows the results of the numerical treatment. The dotted 
curve shows the analytic estimate derived above, where we set 6 = 1.175 to take into account the 
(small) acceleration of the fluid layer, relative to the simple assumptions of the derivation. With 
this correction, the analytic estimate is in good agreement with both the numerical and self-similar 
results. 

These analytic derivations allow us to make an approximate transformation between the effects 
of increasing the initial density [A > 2) and increasing the starting inward velocity (foo > 0). If 
we equate the derived expressions for mo for the two cases (eqs. [56] and [61]), we find that to a 
reasonable approximation A and Voo are related through the cubic equation 

A^{A-l)^A{l + v^f . (62) 

For a typically observed value of v^o = 0.5, for example, the corresponding value of A that produces 
the same infall rate is .4 ~ 2.5. For this case, mo ~ 2.1 and the infall rate is larger than that of 
the infall-collapse solution (Shu 1977) by 115 percent (a factor of 2.15); this result agrees with a 
recent calculation for the collapse of a magnetized singular isothermal toroid (Allen et al. 2003b). 
Observations indicate that inward motions span a range 0.04 - 0.1 km/s (Lee et al. 1998), which 
correspond to either starting speeds in the range Voo = Uin/ag = 0.1 — 0.5 or overdensities in the 
range ^ = 2.1 — 2.5. The largest value of Voo that one should consider is Voo = 1, which limits the 
overdensity parameter to the range A < 2.9 (using eq. [62]). 

4.4. Relation to Observations 

Given the three parameter family of collapse solutions found thus far, we need to specify the 
portion of the parameter space (A, foo 7) that corresponds to observed molecular cloud cores. Since 
the dynamic equation of state (specified by 7) does not affect the mass infall rates as much as the 
other variables (Figures 1 - 7), we focus this discussion on cases with 7 = 1 = F. This paper is 
motivated by the realization that some starless cores are observed to have subsonic inward motions 
(Lee et al. 1999, 2001) at roughly half the sound speed, i.e., nj„ as/2. As mentioned previously, 
these observed motions could arise in two different ways: (1) The inward speed could be part of 
the initial conditions; in this case, we model the collapse using a starting velocity v^o 7^ 0. (2) 
The observed speeds could arise because the core has already begun to collapse on the outside; in 
this case, we model the collapse using overdense initial states (^ > 2), which have zero velocity at 
i = but soon develop subsonic speeds in the outer regions (beyond the expansion wave). We show 
here that both of these possibilities can be made consistent with current observations, although the 
^00 7^ case is favored. 
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Both scenarios (.4 > 2 or > 0) are constrained by the observation that the cores in question 
exhibit inward motions but do not contain detectable infrared sources. The theory is constrained 
because the observed inward velocity must be large enough and the predicted protostellar luminsoity 
must be small enough. For overdense initial states, the predicted inward speed Uin at a given radius 
robs of observation is given by 

Uin = {A- 2)alt/robs , (63) 

where t is the time since the onset of collapse. The luminosity Lth of the central source can be 
written in the form 

GMM ^ 2 

Lth = ^—fT- =^rnl-^t, (64) 

where is the fraction of the total available luminosity Lq (see eq. [1]) that is realized in the 
central source. To match the observations, the predicted inward speeds must be as large as those 
seen, Uin > Uobs- At the same time, the predicted luminosity must be less than the observed limit, 
Lth < Llim (where Lum = 0.1-1 Lq). The first of these constraints implies a lower limit on the 
time t since collapse began, while the second constraint implies an upper limit on t. In order to 
satisfy both constraints, the following bound must be met 

J^—<Liim ^*^ '^ ^ ~0.05 -^ ^ L, (65) 

as Tohsai rn^ rn^ 

where L = Lh^/^ILq), and the numerical constant on the right hand side is obtained by assuming 
robs = 0.1 pc, as = 0.2 km/s, and = 2 x 10^^ cm (Stabler ct al. 1981). The function F{A) = 
{A — 2)/mQ, considered as a function of A, has a maximum value of about 0.112 at v4 ~ 2.5 (found 
by evaluation), and we obtain the constraint 

jT— < 0.006 L. (66) 
ag 

The observed inward velocities lie in the range Uobs/o-s = 0-25 — 0.5, so this constraint requires 
the luminosity to be a small fraction of the total available luminosity, namely T < 0.02. This tight 
constraint greatly limits the available parameter space for the case in which the observed inward 
velocities are part of the collapse (.4 > 2). 

For the scenario in which the inward velocities are part of the initial conditions, the observation 
that Voo = Uin/ as = 0.25 — 0.5 is built into the boundary conditions and does not further constrain 
the solution. However, the requirement that the core does not display a detectable infrared source 
implies a time constraint, 

:Ft < Liin,^ ^ 1.3 X lOVr {L/mD , (67) 

where we have used the same numerical values as before. In order to account for Uobs/O's = 0.25 
- 0.5 = ■Uoo) the implied values of mo lie in the range niQ = 1.6 - 2.2. As a result, we obtain a 
constraint of the form J^t < 5100 yr. In order to allow for a respectably long viewing time, say 
t ~ 5 X 10^ yr (about half the expected collapse time), we would need 0.10. 
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Although both scenarios can be made consistent with observations, the case of overdense initial 
states is much more constrained. This overdense scenario - in which the observed inward motions 
arc part of the collapse itself - requires extremely inefficient power generation in the central source, 
< 0.02. These solutions are also constrained in that they cannot be extended to i < to become 
"complete" solutions (again, see Hunter 1977). In contrast, the Uqo / scenario - in which the 
observed inward motions are part of the initial conditions - is less constrained, but still requires 
.F<0.10 under reasonable assumptions. 



5. COLLAPSE SOLUTIONS FROM NON-ISOTHERMAL INITIAL STATES 

In this section, we consider the collapse of molecular cloud cores with initial density profiles 
that are shallower than p oc (the profile of the singular isothermal sphere). As discussed in 
previous papers (e.g., Lizano & Shu 1989, McLaughlin &: Pudritz 1997), these shallower density 
profiles arise from hydrostatic equilibrium models of gaseous spheres with equations of state that 
are softer than isothermal. In the present context, such states correspond to static F < 1. The self- 
similar collapse of such spheres have been considered previously for the case of 7 = F (see Cheng 
1977, 1978, and the appendices to McLaughlin & Pudritz 1997), although the focus has been on 
hydrostatic initial conditions with v^o = 0. In this section, we generalize the self-similar collapse 
problem to include nonisothcrmal states that arc cither overdense (A > 1) or contain nonzero 
inward velocities (voo > 0). Notice also that the motivation for considering equations of state that 
are softer than isothermal applies primarily to the static equation of state, and not necessarily 
to the dynamic equation of state. We thus generalize these collapse solutions to include the case 

7 7^r. 



5.1. Self-similar Collapse Solutions with a Single Equation of State 

We first consider the case where 7 = F and generalize previous work to include non-zero 
velocities. As discussed earlier, the observed inward motions can arise in two conceptually different 
ways: Such velocities can be "part of the collapse" and hence would imply overdense initial states 
with A > 1; the velocities could also arise as "part of the initial conditions" and would imply A = 1 
and Voo > 0. In either case, we can obtain self-similar solutions using standard methods. The 
equation of motion reduce to the form 

dv y da , , . , , 

where we have defined y = (F — 2)x + v. In order to consider overdense states and nonzero initial 
inward motions, we must specify the outer boundary conditions. Keep in mind that the limit x — ^ 00 
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corresponds to either r ^ oo or t ^ because of the nature of the similarity transformation. The 
starting states are thus a generaUzation of the results given in §4.1 and take the form 



a 



APox-'/(^-^Al - Ao^^^ly x-V(2-r) I , (70) 



^ = -^00 - Aox-r/P-r) , (71) 
where we have defined the dimensionlcss coefficients 



2r(4 - 3r) 
(2 - r)2 



i/(2-r) 



and Ao = A/3o [l - A"^'-^)] . (72) 



Notice that inclusion of the constant term Voo in the velocity expansion breaks the self-similarity. 
Solutions with Voo 7^ can be found, but, unlike the case of isothermal solutions discussed in the 
previous section, the results depend on the choice of outer boundary condition (sec below). 

Wc note that the equations of motion can contain critical points, which can be considered as 
generalized sonic points. The condition for the flow to pass through a critical point take the form 

[(r - 2)x + vf = Fa^-i . (73) 

Although the fluid flelds are continuous at the critical points (deflned by eq. [73]), the derivatives 
(here, with respect to the similarity variable x) of the fluid fields need not be. As a result, care 
must be taken when numerically integrating through critical points. One can expand the functions 
in the neighborhood of the critical points and use the resulting analytic forms to integrate through 
the point. This treatment has been discussed previously in the literature (e.g., see Shu 1977 and 
McLaughlin &; Pudritz 1997) and need not be belabored here. In the present application, only 
solutions starting from hydrostatic initial conditions actually go through the critical points. For 
starting states that are even slightly overdense, or ones that start with small initial inward velocities, 
the solutions avoid the critical points altogether. 

The basic results of this section are shown in Figure 8 and 9. Keep in mind that the infall rate, 
the quantity of interest in this collapse calculation, is given by the analytic expression of equation 
(31), with the reduced mass rriQ determined numerically (as shown in the Figures). Figure 8 shows 
the reduced mass mo as a function of the index F = 7 for varying values of the overdensity parameter 
A. The results are plotted as the ratio mo/mH, where mn is the mass enclosed within the expansion 
wave for the case of hydrostatic equilibrium, i.e., with A = 1. During collapse, the overdense states 
develop a rapid inward flow well beyond the nominal location of the expansion wave (which is most 
meaningful for the hydrostatic case). As a result, the mass infall rates become rather large, and 
vary with overdensity more sensitively than in the case of the isothermal collapse calculations of the 
previous section. Specifically, as the density increases by a factor of 2 (A = 1 ^ 2), mo increases by 
a factor of ~1840 in the logatropic limit F — > 0, compared to only a factor of 5.7 in the isothermal 
limit F — >^ 1. Notice also that because of the rapid speed at large r, the ratio mo/mH can be greater 
than unity. 
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Figure 9 shows the corresponding results for the reduced mass mo for varying values of the 
fixed inward speed foo- As before, the results are presented in terms of the ratio mo/m//, where uih 
is the hydrostatic value (which is constant for all t^oo)- These solutions have a more complicated 
interpretation than their isothermal counterparts of the previous section. For isothermal initial 
conditions, the reduced speed v is just the physical speed u scaled by the (constant) isothermal sound 
speed. For nonisothermal conditions, the physical speed u = wv{x), where v is the (dimensionless) 
reduced velocity and w = k^^'^ (iTrGt'^Y^^^^^'^ is (esssentially) the effective sound speed at the head 
of the expansion wave. Unlike the isothermal case, however, the scaling speed w is time dependent. 
As a result, the solutions arc not self-similar and depend on the boundary conditions, i.e., on the 
starting value xq of the similarity variable where integration of the equations of motion begins. 
Nonetheless, solutions can still be found (see Figure 9), but the interpretation is not as clean as 
before. The solutions shown in Figure 9 have been calculated by applying the outer boundary 
condition (eq. [71]) at xq = 100. In contrast to the self-similar isothermal case of the previous 
section, the resulting values of tuq depend on the choice of xq. Nonetheless, the curves shown in 
Figure 9 ilhistratc the expected property that collapse solutions with initial inward motions lead 
to larger infall rates than those found in the hydrostatic limit. 



5.2. Estimate for mo for Hydrostatic Initial Conditions 

In this section, we derive an estimate for the reduced mass rriQ at the origin for the case of 
hydrostatic initial states. This estimate is found by calculating the value of mo that would arise 
for a pressure-free collapse. Since the pressure forces are small, this calculation should provide a 
reasonable estimate of the true value of mo- The pressure forces are not zero, however, and retard 
the flow so that the true value of mo will be smaller than given by this estimate. We correct for 
this difference using a simple analytical fit. 

To start, we consider a fluid layer at an initial radius tq. Because of the self-similar nature of 
the problem, the value of vq does not matter. For the collapse solutions of interest, the fluid layer 
will remain at rest until the expansion wave reaches the radius ro; afterward, the layer begins to 
fall. As the layer falls inward, the total mass enclosed within its current radial position remains 
constant and is given by the hydrostatic profile (§3.2) evaluated at ro, i.e., M(ro). 

The time tf required for a fluid layer to fall from an initial radius ro is given by 

tf = ^rl/\2GM)-'/\ (74) 

where we have assumed no pressure forces so that this time scale is a lower limit. Notice that if 
the fluid layer is not assumed to start at rest, but rather has an initial inward velocity vq, then the 
leading numerical coefficient will be different (not tt/2). 



The time tc for the expansion wave to reach the initial radius ro is given by 

tc=[Aro/xHf^''-^\ 



(75) 
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where A is the coefficient in the similarity transformation given by equation (22). When the fluid 
layer reaches the origin, the total time tr that has elapsed is given by the sum 

tT = tc + tf. (76) 

The total mass in the central object (at the origin) is the mass enclosed within the radius tq. 
The total mass enclosed with the expansion wave is given by the hydrostatic mass profile (see §3.2) 
evaluated at the location of the expansion wave radius at time tx, i.e., at 

ru = xhA-H^-^ . (77) 



Thus, the ratio of masses is given by 

(4-3r)/(2-r) 



ruH 



rH 



to 

tT 



(4-3r) 

(78) 



In this case, since the initial states are hydrostatic, tq = ru and the ratio mo/mH is always less 
than unity (compare with the overdense states of §4.3). 

We can use the similarity transformation and the hydrostatic initial state to write the pressure- 
free collapse time tp as a fraction of the total time tx- We also introduce a correction factor g > 1, 
which can be a function of the equation of state, and which accounts for the fact that the true 
infall time is longer than the pressureless free fall time calculated above. After some algebra, the 
mass ratio can be written in the form 

( >! -(4-3r) 

^= i + 5l(7/r)^/^(2-r)-v^ . (79) 

Since we already have a closed form expression for the enclosed mass niH (equation [30]), we thus 
have a completely analytic estimate for mo. Also, since this estimate was derived assuming no 
pressure forces, this expression (with g = I) represents an upper limit for the reduced mass mo. As 

an example, we consider the limit of the singular isothermal sphere (Shu 1977) where 7 = F = 1. 
For this case, rriH = 2 and wc obtain the estimate mo = 8/ (4+7r) « 1.12, which compares reasonably 
well with the calculated value of mo = 0.975. Although the expression (79) is only accurate at the 
10 percent level, it illustrates how the mass infall rate (which is proportional to mo - equation [31]) 
varies with the equation of state (mostly the static index F). 

Figure 10 shows the values of the central mass mo (which sets the mass infall rate) as a function 
of the index F in the static equation of state for the simplest case in which the dynamic 7 is the 
same as the static F. The solid curve shows the results obtained from solving the self-similar 
equations of motion (this result has been obtained previously - see McLaughlin & Pudritz 1997). 
The dashed curve shows our analytic approximation for mo, where we take into account the fact 
that pressure slows down the infall so that the time required to fall from a given starting radius 
is longer by a factor of g (see equation [79] ) . This effect is greater for softer equations of state, so 
we take g = g{T) = (4/3) (2 — T)^/^^, where the index of 9/10 is found by fitting the numerically 
determined curve. The result provides almost an exact fit. 
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5.3. Different Dynamic and Static Equations of State 

Finally, we consider collapse flow with nonisothermal initial states and varying choices for the 
dynamic equation of state. Figure 11 shows the result of keeping a fixed static T while allowing the 
dynamic equation of state (dynamic 7) to vary. As shown, the variation in the dynamic equation 
of state has relatively little effect on the resulting reduced masses uiq, except when the dynamic 
index 7 is much larger than the static index T. In this latter regime, the pressure forces within the 
incoming flow become large and the infall rates arc suppressed. This suppression is best illustrated 
by considering a standard reference case where the static index F is arbitrary, but the dynamic 
equation of state is chosen to be isothermal and fixed (7 = 1). Figure 12 shows the result for 
varying static indices F and a fixed (isothermal) dynamic equation of state. The ratio mo/mH of 
the central mass variable to the mass within the fiducial expansion wave (this ratio determines the 
mass infall rate) is plotted as a function of the static index F (as before, the denominator uses the 
hydrostatic value of rriH with A = 1). The solid curves show the case of 7 = F for two choices of 
the overdensity (A = 1.5 and 2). The dashed curves show the resulting values of mo/mH for the 
same starting conditions, and the same static F, but with dynamic 7 = 1. This generalization to 
7 7^ F leads to increased pressure forces during the collapse (relative to the case of F = 7) and the 
resulting values of mo are smaller. The discrepancy grows with decreasing static index F. 

Before leaving this section, we note that the mathematical formulation of the collapse problem 
used in this paper (§2) docs not allow self-similar solutions for the case with F ^ and 7 7^ F. 
In other words, in the logatropic limit, collapse with other dynamic equations of state remains 
undefined. One manifestation of this problem can be found in equation (16) for the reduced 
pressure. In the limit F — > (7 7^ 0), j3 ^ constant instead of the required limiting form p = log a. 
The class of solutions found in this paper are known mathematically as self-similar solutions of the 
first kind (Barenblatt 2003). Although the case of F 0, 7 7^ does not allow for a self-similar 
solution of the first kind, it remains possible that the problem will admit self-similar solutions of the 
second kind, also known as incomplete similarity (Barenblatt 2003; Zeldovich 1956; von Weizsacker 
1954). This complication is beyond the scope of this paper and will be left for future consideration. 

6. CONCLUSION 

6.1. Summary of Results 

In this paper, we have explored a wide variety of collapse solutions for the collapse problem 
in star formation. These solutions describe the collapse of idealized molecular cloud cores and/or 
the collapse of kernels, the subcondensations within large cores that form stellar groups and clus- 
ters. Recent observations of star forming regions (e.g., Lee et al. 1999, 2001) indicate that the 
outer boundary condition for the collapse flow is not static; instead, the collapse flow must match 
smoothly onto a subsonic, but finite, velocity field at large radii. This paper shows that a wide 
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variety of self-similar solutions have small but finite velocity fields at large radii, and we argue that 
solutions of this type provide a better description of protostellar collapse. The most important 
ramification of this change in the outer boundary condition is that the mass infall rates will be 
larger for cores with nonzero flows. We have quantified this change using the semi-analytic tech- 
nique of finding similarity solutions; we have verified this result using a full numerical treatment 
of collapse and have derived analytic bounds/estimates to provide greater physical understanding 
of the process. All three methods are in good agreement and show how the infall rates depends on 
the parameters of the problem. Our specific results are listed below: 

[1] In general, the equation of state that describes the thermodynamics of the collapse flow - that 
which determines the evolution of entropy - need not be the same as the effective equation of 
state that sets up the initial conflguration. For the case of polytropic equations of state, this 
generalization means that the dynamic index 7 need not be the same as the static index T. This 
paper shows that similarity solutions can be found for the more general problem containing two 
equations of state. In addition, although previous work has focused on collapse flows that match 
onto static outer boundary conditions, this paper shows that gaseous spheres with inward flows at 
large radii also can be described by similarity solutions; we parameterize this latter generalization 
by introducing the reduced speed v^c (the speed at a; ^ 00). Including the ovcrdensity parameter 
A, which has been considered previously, this paper shows that self-similar collapse solutions exist 
for a four parameter family of conditions, where points in this parameter space can be written in 
the form (F, 7, A, Uoo). 

[2] The initial hydrostatic equilibrium configuration is determined by the static equation of state. 
In order for equilibrium configurations to exist, the index of the static equation of state must satisfy 
the constraint F < 4/3. The power-law index jj, of the density profile is given by 

A = -^. (80) 

This result shows that the equilibrium density distribution becomes less centrally concentrated as 
the static equation of state becomes softer (as the static index F decreases). 

[3] The time dependence of the mass infall rate also depends only on the static equation of state 
and is given by 

M ~ t3(i-r) _ (g;^) 

As a result, static equations of state that are stiffer than isothermal produce mass infall rates that 
decrease with time. Similarly, static equations of state that are softer than isothermal lead to 
mass infall rates that increase with time. Since we expect the equilibrium configurations of actual 
cloud cores to correspond to the softer equations of state, the mass infall rate should either remain 
constant (as for the isothermal case) or increase with time. This result holds for any value of 
dynamic 7. 

[4] Collapse solutions which approach a ballistic free-fall form in the inner region {x ^ 1) exist for 
all dynamic equations of state with indices 7 < 5/3. For these solutions, the reduced fluid flelds 
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approach the asymptotic forms 

m~mo, u~-V2^x"^/^ a ~ (4 - Sr) v^m^a;"^/^ (82) 

in the hmit x ^ 0. The reduced mass mo sets the magnitude of the mass infall rate. As shown 
by equation (31), the infall rate is determined, except for the reduced mass, by the similarity 
transformation. 

[5] For collapse starting from an isothermal static equation of state (which implies p ~ r~^), we 
have solved the full equations of motions (the original partial differential equations) in addition to 
the similarity equations. These numerical results provide important confirmations of the similarity 
approach: (A) The full numerical treatment verifies the solutions found semi-analytically; in par- 
ticular, they predict the same mass infall rates. (B) The numerical results show that non-singular 
initial configurations, in particular those that have a central region of constant density, approach 
the similarity solutions asymptotically in time. Moreover, the convergence occurs relatively rapidly 
and smoothly, especially for the obscrvationally preferred cases with A > 1 and Vr^ ^ 0. 

[6] We have derived analytic estimates for the reduced mass mo which sets the mass infall rate 
through equation (31). For isothermal starting states, we obtain analytic expressions for tuq as a 
function of overdensity A (see equation [56] ) and as a function of the initial inward speed f oo (see 
equation [61]). These expressions predict the same infall rates calculated from both the numerical 
treatment (solving the partial differential equations) and our semi-analytic approach (the similarity 
solutions) as shown in Figures 6 and 7. For the case of non-isothermal initial conditions, we have 
derived an analytic expression for mo as a function of the index F (sec equation [79]) for collapse 
with 7 = F. The resulting expression is in good agreement with the values calculated from the 
similarity solutions (Figure 10). 

[7] For collapse solutions with 7 7^ F, the infall rates are esssentially the same when the two indices 
are comparable (see Figures 1 and 11). For soft equations of state F <C 1, however, the infall rates 
are greatly suppressed for relatively large values of dynamic 7 ~ 1 (see Figure 12). In the logatropic 
limit for the static equation of state, F — >^ 0, this formulation does not admit self-similar solutions 
of the first kind for arbitrary dynamic 7. 

[8] The observational motivation for this study is that starless cores can exhibit inward motions at 
about half the sound speed, nj„ ^ as/2. These motions could arise from two conceptually different 
sources. If the core is overdcnse A > 1, then the entire core - including the outer regions - will begin 
to collapse, with observable inward motions appearing on the outside before an observable infrared 
source appears in the center. Alternatively, the inward motions could be part of the initial conditions 
- they could be present at the effective zero point of time (the start of collapse) , presumably left over 
from the core formation process. Both cases are compatible with current observations, although 
the overdense scenario is more constrained. The net result, for both scenarios, is that the infall 
rate with nonzero initial velocities is larger (than in the hydrostatic case) by a factor of two (see 
also Allen et al. 2003b). 
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6.2. Discussion 

The specific results outlined in §6.1 represent a significant mathematical generalization of the 
self-similar collapse problem, which now includes the four parameter family of solutions given by 
(r,7,A, Voo). Moreover, the analytic expressions derived for mo(A), mo('yoo)) and mo(r) provide 
both a convenient means of obtaining the infall rates (through equation [31]) and a path toward 
greater physical understanding of the infall process. This work has additional astronomical impli- 
cations, however, as we discuss next. 

The generalized collapse solutions presented here are motivated by astronomical observations 
which indicate that cores have inward motions at large radii at the start of the collapse phase (Lee 
et al. 1999, 2001). The first implication of these new solutions is that molecular cloud cores will 
collapse with larger mass infall rates due to these inward velocities. When matched to the observed 
boundary conditions, our collapse solutions imply that the infall rates will be a factor of two larger 
than previous estimates. 

Even with the previous (smaller) infall rates, observed protostellar objects have smaller lumi- 
nosities than those predicted by theoretical infall rates (e.g., Kenyon & Hartmann 1995), unless a 
great deal of the energy is stored as orbital motions. With the higher infall rates now indicated, 
this luminosity problem is worse by a factor of two, which means that even more energy must be 
stored as rotational energy in the disk. With increased infall rates, the predicted time scales for 
star formation are shorter by this same factor of two. As argued from both observational (Myers 
k, Fuller 1992) and theoretical (Adams & Fatuzzo 1996) considerations, stars form on a time scale 
tsf ~ 10^ years over a wide range of stellar masses. Although the mean time is shorter, the pre- 
diction that the distribution of time scales is much narrower than the distribution of stellar masses 
(the IMF) remains the same. 

The theoretical picture of star formation developed in the 1980s considered molecular cloud 
cores to be (nearly) hydrostatic at the beginning of the collapse phase. More recent observations 
of the numbers of starless cores (e.g., Jijina et al. 1999) indicate that core formation must take 
place more rapidly and is unlikely to involve a purely quasi-static process. Some modification to 
the picture is thus necessary. Many recent authors (see the review of MacLow & Klessen 2004) 
have been exploring the opposite limit in which the star formation process is fully dynamic and 
that cloud cores never really form at all (at least not as separate, physically well-defined entities). 
The picture of star formation emerging from this study is intermediate between these two extremes: 
Cores form relatively quickly (through physical processes not calculated here), but still represent 
the initial conditions for the following stage of protostellar collapse. With their rapid formation, 
the cores display non-zero inward motions at large radii as observed. When this dynamic element 
is incorporated into the collapse solutions (as calculated in this paper), the subsequent collapse 
proceeds much as in previous models but with a larger infall rate (due to the head start provided 
by the initial velocities). As a result, the scenario retains the successes of the standard picture 
(Shu et al. 1987), such as the spectral energy distributions of young stellar objects (Adams et al. 



-29- 



1987) and their corresponding observed emission maps (Walker et al. 1990). At the same time, the 
picture can be made consistent with faster formation times for cfoud cores (Jijina et al. 1999) and 
the observed inward motions at large radii (Lee et al. 1999). 

This modification to the initial conditions - inward velocities due to either overdense starting 
states or nonzero as an initial condition - alleviates one theoretical difficulty associated with 
protostellar collapse, but raises another. On the positive side, the inward velocities ensure that 
non-singular cores collapse smoothly, and more rapidly approach the self-similar solutions. On the 
negative side, given their dynamic formation, the cores no longer display a clean separation between 
their formation stage and their collapse stage. 

In the (old) hydrostatic picture, molecular cloud cores live on the edge between collapse and 
expansion. For example, in the isothermal case (Shu 1977) where po = Aa'^/4:TrGr^ and A = 2 
corresponds to equilibrium, cores with A> 2 will collapse but cores with A < 2 will actually expand! 
If cores were formed in states indistinguishable from the hydrostatic equilibrium configuration with 
A = 2, then they would always be in danger of expanding. In retrospect, it seems unlikely that 
the fate of a cloud core - to collapse or to expand - should depend with hair-trigger sensitivity on 
whether ^ = 2 + eor^ = 2 — e (where e ^ 1). In contrast, cores with initial inward velocities (as 
now observed) will collapse under a robust set of circumstances. 

The initial inward velocities also affect the issue of cores with constant density centers (as 
observed by Ward-Thompson et al. 1999, Bacmann et al. 2000, Tafalla et al. 2004). For fiat-topped 
cores that start close to hydrostatic equilibrium, realistic collapse flows approach the self-similar 
form when the ratio rout/fc ^20 (Foster & Chevalier; Whitwortli et al. 1996). Roughly speaking, 
the central (constant density) region collapses first, on a time scale given by At « (37r/32Gpo)^/^ 
(e.g., Spitzer 1978). For later times t > At, the flow approaches the self-similar form as long 
as more of the cloud exists to subsequently collapse (this condition holds for rout/f'c^2Q). For 
cores with initial velocities, however, the time required for the central (constant density) region to 
collapse is reduced (as calculated in equation [52]) so that the requirement on the ratio rout/Tc is 
less constraining. In order words, because observed cores have initial inward velocities, the self- 
similar solutions provide a better description of the collapse flow (than if the cores were exactly 
hydrostatic) . 

An explicit demonstration of this trend is shown in Figure 13, which plots the mass infall rates 
for a set of flat-topped cores as a function of time. This set of numerical collapse calculations begins 
with isothermal cores with sound speed ag = 0.2 km/s and an extended central region of nearly 
constant density spanning rc ~ 10^^ cm. The initial states are slightly overdense, with A = 2.2. 
The solid curve shows the resulting mass infall rate (expressed in Mq yr~"^) for a starting state 
with no initial velocity. As found previously, the infall rate is essentially zero to begin with, and 
then exhibits a sharp spike at the time when the constant density region reaches the origin; the 
infall rate then decreases with time and eventually approaches the asymptotic value M = niQa^/G 
predicted by the similarity solution. The dotted (dashed) curve shows the corresponding mass infall 
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rates for starting inward speeds of Uin = 0.5 as {uin = as)- When Uin ^ 0, the mass infall rates 
reach their peak values sooner and they more rapidly approach their asymptotic values (which arc 
within about 10 percent of those predicted by the self-similar calculations). In addition, for cores 
with Uin 7^ 0, the maximum infall rates are smaller and the asymptotic values are larger, so the 
time variability is less severe. More specifically, the ratio of the peak infall rate to its asymptotic 
value is ~8.4 for nj„ = and only 2 - 3 for Uin ^ 0. It has been suggested (e.g., Ward-Thompson et 
al. 1999) that the spike in the infall rate corresponds to the Class O phase of protostellar evolution. 
Although this work suggests that initial inward velocities smooth out the spike to some degree, the 
infall rates are nonetheless larger early on and can correspond to the Class phase. 

The new difficulty raised by the presence of initial inward speeds has two components. First, the 
physical processes that explains core formation must be able to account for the observed speeds. 
Second, the core formation epoch must be more dynamic than considered previously, and the 
distinction between the core formation stage of evolution and the core collapse stage is not as 
clean. While this paper clarifies the collapse phase of star formation, the (earlier) core formation 
phase requires further theoretical work. 

Finally, one must keep in mind that this work generalizes spherical collapse solutions in the 
absence of magnetic fields and angular momentum. Nonrotating solutions are expected to apply 
at larger radii where the flow is nearly spherical. The effects of rotation become important in the 
inner part of the flow, where our solutions can be matched onto inner solutions as outlined in §3.5. 
The inclusion of magnetic fields into the collapse problem has been studied in a variety of related 
contexts (see GaUi & Shu 1993ab; Li & Shu 1996, 1997; Allen et al. 2003ab; Shu et al. 2004), 
including one collapse calculation with nonzero initial velocities (see Figure 7 of Allen et al. 2003b). 
With the inclusion of rotation (at rate $7) and magnetic fields B, one is left with a six parameter 
family of initial conditions (F, 7, A, Wqo, f^, B), where the velocity field Voo(r), the rotation profile 
r2(r), and the magnetic field strength B(r) can all be functions of the spatial coordinates. With 
this paper (which explores collapse with Voo / and 7 / F), the effects of these parameters on 
collapse have all been considered individually. The challenge for the future is thus to identify the 
portion of this parameter space that applies to protostellar collapse and to properly incorporate all 
of the relevant effects (simultaneously) into the collapse solutions. 
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APPENDIX: EXISTENCE OF FREE-FALL SOLUTIONS 

In this appendix, we show that baUistic free-fall solutions exist in the inner regime of the 
collapse flow for any dynamic equation of state with index 7 < 5/3. 

We begin by taking the limiting form of the equations of motion for the case 

x^\v\, x<l, (83) 

which defines the inner regime of the collapse flow. In this limit, the general expressions for the 
mass m{x) and the pressure p{x) reduce to 

, , vx'^a . 

m(x) = , (84) 

^ ^ 3a + 2 ' ^ ' 

and 

p{x) = Coa' {avx^/ay = da^ . (85) 
Next, we note that the final equation of motion takes the form 

a ax V ax x 
which can be immediately integrated to obtain 

avx'^ = constant . (87) 

Comparing this result with the limiting form for the mass m(x) shows that the enclosed mass 
always approaches a constant value in the inner limiting regime, i.e., 

■m{x) — mo = constant . (88) 

This result is quite general; it was obtained by assuming only that the coordinate x is small (which 
defines what we mean by the inner regime of the flow) and by assuming that x <C \v\. This latter 
assumption implies that we are considering only infall solutions in this argument. As a general rule, 
outflow solutions to the equations of motion will also often exist. The outflow solutions generally 
have V ^ and x ^ and hence the results of this current argument are not applicable. 

Given that the reduced mass approaches a constant value in the inner regime, the reduced 
pressure can be simplifled also, i.e., 

p{x) Ca^ , (89) 

where we have deflned a new constant C. This result makes sense: it means that the pressure is 
given by a pure polytropic form in the inner limit, where the polytropic index is determined by 
the dynamic equation of state. In this limit, the gas has lost all memory of the original "static" 
equation of state. 
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The equation of motion (the force equation) reduces to the form 



i;— + 7CaT-2— + ^ = 0. 
ax ax 



Our objective is to show that a self-consistent free-fall solution exists for certain values of the index 
7. If a free-fall solution exists, then the velocity must approach the form 

v = -{2mo)^/^x-^/\ (91) 

where the minus sign indicates that the material is falling inward. Using this form for the velocity 
in the equation of motion (90), wc sec that the first velocity term exactly balances the third 
gravity term. This solution will exist provided that the second pressure term vanishes. Since 
m —>■ mo ^ vx^a in this limit, a ~ x^^l"^. The ratio Tl of the pressure term to the gravity term 
thus becomes 

7^ ~ constant x x^^'^^^l'^ . (92) 
This ratio vanishes in the limit of small x provided that 

7 < 5/3 . (93) 

Thus, self-consistent free-fall solutions exist in the inner limit for all dynamic equations of state 
with indices 7 < 5/3. 



NOTE ADDED IN PROOF 



Related work by McKee &; Holliman (1999) shows that singular isothermal spheres (with static 
r = 1) that are initially in hydrostatic equilibirum will not collapse if dynamic 7 > 4/3. On the 
other hand, this paper finds self-similar collapse solutions for 7 < 5/3 (e.g., see Figure 1). The 
solutions presented here use starting states that are overdense (not in hydrostatic equilibrium) so 
they will always begin to collapse. In order for continued collapse to take place, the gas must 
be able to get to the inner, free-fall regime before pressure forces take over. In the usual star 
formation scenario, the gas strikes a stellar surface, releases energy, and collapse can continue. In 
the physical (not self-similar) problem, inner boundary conditions thus help determine whether 
or not sustained collapse takes place. As a result, for cases with dynamic indices in the range 
4/3 < 7 < 5/3, the self-similar collapse solutions of this paper may not always provide good models 
for physical collapse solutions that include additional scales (e.g., finite mass, fixed outer boundary, 
finite central density). The determination of which physical conditions allow continued collapse 
in this regime must be left for a future research project. However, the majority of our collapse 
solutions, those with dynamic 7 < 4/3, are not affected by this issue (we thank Chris McKee for 
useful discussions on this topic). 
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Fig. 1. — For initial states with isothermal statie F = 1, the mass infall rate M = moa^/G. This 
plot shows the variation of mo with the overdcnsity of the initial configuration for various values 
of the dynamic 7. Hydrostatic equilibrium corresponds to „4 = 2. The upper dashed curves show 
results for 7 = 0.6 and 0.8 (7 < 1); the lower dotted curves show results for 7 = 1.2, 1.4, and 1.6 
(7 > 1); the solid curve corresponds to 7 = 1 (see Shu 1977). 
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Fig. 2. — For initial states with isothermal static F = 1, the mass infall rate M = moa^/G. This 
plot shows the variation of ttiq with the initial inward velocity Voo = Uin/ag for various values of 
the dynamic 7. Hydrostatic equilibrium corresponds to Uin = 0. The upper dashed curves show 
results for 7 = 0.6 and 0.8 (7 < 1); the lower dotted curve shows results for 7 = 1.2; the solid curve 
corresponds to 7 = 1. 
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Fig. 3. — This figure shows the numerically determined collapse solution for a gaseous sphere 
with the density appropriate for hydrostatic equilibrium {A = 2) and an initial inward velocity of 
magnitude Uin = O.bag. The top panel shows the resulting density profile at four times (scaled to 
typical parameters for a molecular cloud core). The bottom panel shows the corresponding inward 
velocity profile at the same four epochs. Notice that v — ^ constant in the limit r — ^ oo and the 
profiles have the same (self-similar) form at smaller radii. 
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Fig. 4. — This figure illustrates how the numerically determined collapse solution for a gaseous 
sphere approaches a self-similar form. The numerically determined solutions for both the density 
(top panel) and velocity (bottom panel) are rescaled according to the (expected) similarity trans- 
formation and then plotted on top of each other. The result is one smooth function (for each panel) 
and hence the numerically determined solutions at any given time are a stretched version of the 
solution at an earlier time, in keeping with the expectation of self-similarity. The corresponding 
self-similar solution, shown by the filled circles, is in good agreement. 
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Fig. 5. — Numerically determined infall rates shown as a function of the overdensity parameter A 
and the initial (inward) velocity. All collapse models use isothermal conditions so that 7 = 1 = T. 
Each curve shows the numerically calculated values of niQ as a function of A for a given starting 
speed Uin- The lower curve shows the result for Uin = 0. The next four curves, in ascending order, 
show the results for Uin/ag = 0.25, 0.5, 0.75, and 1.0, where Ug is the isothermal sound speed. 
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Fig. 6. — Comparison of self-similar, analytic, and numerical determinations of mo, the dimension- 
less mass that specifies the infall rate for isothermal initial conditions. The analytic estimate for 
tuq is shown by the dotted curve; the values calculated from the self-similar equations of motion are 
shown by the solid curve; the results from a numerical treatment (solving the partial differential 
equations) are shown as the dashed curve. All models have initial states with isothermal static T 
= 1 and dynamic 7 = 1. Keep in mind that for all of these configurations, the mass infall rate is 
given by M = moa^/G so that mo specifies the infall rate. 
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Fig. 7. — Comparison of self-similar, analytic, and numerical determinations of tuq for varying 
initial inward velocities Uin and for isothermal initial conditions. The mo values calculated from 
the self-similar equations of motion are shown by the solid curve; the results from a numerical 
treatment (solving the partial differential equations) are shown as the dashed curve. The analytic 
estimate for rriQ is shown by the dotted curve; to obtain this form, we assume a fitting parameter 
that takes into account a slight (20 percent) acceleration of the flow before the expansion wave 
passes. All models use initial states with isothermal static F = 1 and dynamic 7 = 1. For all cases, 
the mass infall rate is given by M = moa^ /G so that mo specifies the infall rate. 
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Fig. 8. — Reduced masses uiq as a function of equation of state index T = 7 for varying values 
of the overdensity. The results are expressed as the ratio mo/niH, where uih is the reduced mass 
enclosed within the expansion wave for the case of hydrostatic equilibrium (the reduced masses 
mo for each overdensity is thus compared to the same mn)- The lowest curve shows the case of 
hydrostatic equilbrium with A = 1. The subsequent curves (proceeding upward) show the results 
for A = 1.05, 1.10, 1.25, 1.50, 1.75, and 2.0. Notice that the infall rate, which is proportional to 
mo, varies much more rapidly with the overdensity A for smaller values of index F. 



-44- 




Fig. 9. — Reduced masses mo as a function of equation of state index F = 7 for varying values of 
the inward speed (where v(x) — >■ Voo at large x). The results are expressed as the ratio mo/mH, 
where mn is the reduced mass enclosed within the expansion wave for the case of hydrostatic 
equilibrium (the reduced masses mo for each case is compared to the same mn)- The lowest curve 
shows the case of hydrostatic equilbrium with Vq^ = 0. The subsequent curves (proceeding upward) 
show the results for Voo = 0.1, 0.2, 0.4, 0.6, 0.8, and 1.0. Because the inclusion of v^o 7^ breaks 
the similarity, only the lowest curve is self-similar. The other curves of mo/m//(F) depend on the 
value xo at which the outer boundary condition is applied; xo = 100 for the cases shown here. 
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Fig. 10. — Comparison of the analytic estimates for mo (dashed curve) with the values calculated 
from the self-similar equations of motion (solid curve) for initial states with varying values of the 
polytropic index. In this case, the dynamic 7 is set equal to the static F. The mass infall rate M is 
proportional to mo (although the scalings depend on 7 = F). The dashed curve shows the analytic 
result of equation (79) with ^(r) = (4/3) (2 — F)^/^°; this form takes into account the result, found 
numerically, that polytropes with softer equations of state approach free-fall conditions more slowly 
as F decreases. 
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Fig. 11. — The effect of dynamic 7 on infall rates. The various curves show the value of the reduced 
mass niQ (which sets the infall rate) as a function of dynamic 7 for fixed values of the static equation 
of state, i.e., fixed values of static F (as labeled). The variation in dynamic 7 has relatively little 
effect unless 7 » F; in this latter regime, the reduced masses (and infall rates) become much 
smaller than in the case of 7 = F. 
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Fig. 12. — Reduces mass ratio for varying static indices F and fixed (isothermal) dynamic 7 = 1. 
These curves show the numericahy determined values of mo/niH, which set the mass infall rates, 
as a function of static index F (the denominator uses the hydrostatic value of mn with A = 1). 
The solid curves represent models in which the dynamic index is the same, i.e., 7 = F. The dashed 
curves show cases where dynamic 7 = 1. In other words, the polytropic spheres are allowed to have 
varying initial density profiles (as set by varying F), but the entropy evolution of the gas takes place 
under isothermal conditions 7 = 1. The upper curves correspond to initial states with overdensity 
A = 2; the lower curves use overdensity A = 1.5. 
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Fig. 13. — Mass infall rates for flat-topped cores as a function of time for varying initial infall 
speeds. This set of numerical collapse calculations begins with typical parameters observed in 

flat-topped cores with sound speed as = 0.2 km/s and a central region of nearly constant density 
spanning vq ~ 10^^ cm. The initial states are overdcnse at the 10 percent level, so that A = 2.2. 
The solid curve shows the resulting mass infall rate (expressed in Mq yr~^) for the case with no 
initial velocity. The dotted (dashed) curve shows the corresponding mass infall rate for starting 
inward speeds of Uin = 0.5 as [uin = as)- When Uin 7^ 0, the mass infall rates reach their peak 
values more quickly, the peak values arc smaller, and the infall rates more rapidly approach their 
asymptotic values (as predicted by the similarity solutions). 



